﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using dnAnalytics;
using dnAnalytics.Math;
using System.IO;
using System.Runtime.Serialization.Formatters.Binary;
using System.Drawing;

namespace CalcSharp.Core.Containers
{
    [Serializable]
    public class FloatMatrix : AMatrix
    {
        private readonly float[] mData;
        internal float[] Data { get { return mData; } }

        #region Constructors
        public FloatMatrix(int rows, int columns): base(rows, columns)
        {
            mData = new float[rows * columns];
        }
        public FloatMatrix(int rows, int columns, float value): this(rows, columns)
        {
            for (int i = 0; i < mData.Length; i++)
            {
                mData[i] = value;
            }
        }
        public FloatMatrix(FloatMatrix other): this(other.Rows, other.Columns)
        {
            if (other != null)
            {
                Buffer.BlockCopy(other.mData, 0, mData, 0, mData.Length * sizeof(float));
            }
        }
        public FloatMatrix(float[,] array): this(array.GetLength(0), array.GetLength(1))
        {
            for (int j = 0; j < Columns; j++)
            {
                int jIndex = j * Rows;
                for (int i = 0; i < Rows; i++)
                {
                    mData[jIndex + i] = array[i, j];
                }
            }
        }
        #endregion

        public override void CopyTo(AMatrix target)
        {
            if (target == null) throw new ArgumentNullException("target");
            if (Rows != target.Rows || Columns != target.Columns) throw new Exception("Target and Source row/column count mismatch");

            FloatMatrix Dbl = target as FloatMatrix;
            if (Dbl == null) base.CopyTo(target);
            else Buffer.BlockCopy(mData, 0, Dbl.mData, 0, mData.Length * sizeof(float));
        }

        public override AMatrix Clone()
        {
            return new FloatMatrix(this);
        }

        #region Abstract implementations
        protected internal override void ValueAt(int row, int column, double value)
        {
            float v;
            if (value > float.MaxValue) v = float.MaxValue;
            else if (value < float.MinValue) v = float.MinValue;
            else if (value == double.NaN) v = float.NaN;
            else if (value == double.NegativeInfinity) v = float.NegativeInfinity;
            else if (value == double.PositiveInfinity) v = float.PositiveInfinity;
            else v = (float)value;
            mData[column * Rows + row] = (float)v;
        }

        protected internal override double ValueAt(int row, int column)
        {
            return mData[column * Rows + row];
        }

        protected internal override AMatrix CreateMatrix(int numberOfRows, int numberOfColumns)
        {
            return new FloatMatrix(numberOfRows, numberOfColumns);
        }

        public override double Determinant()
        {
            Compute();
            if (mIsSingular) return 0;
            if (mDeterminant == double.MinValue)
            {
                lock (mMatrix)
                {
                    mDeterminant = 1.0;
                    for (int j = 0; j < mMatrix.Rows; j++)
                    {
                        if (mPivots[j] != j) mDeterminant = -mDeterminant * mMatrix.ValueAt(j, j);
                        else mDeterminant *= mMatrix.ValueAt(j, j);
                    }
                }
            }
            return mDeterminant;
        }

        public override AMatrix Transpose()
        {
            FloatMatrix ret = new FloatMatrix(Columns, Rows);
            for (int j = 0; j < Columns; j++)
            {
                int index = j * Rows;
                for (int i = 0; i < Rows; i++)
                {
                    ret.mData[i * Columns + j] = mData[index + i];
                }
            }
            return ret;
        }

        public override AMatrix Inverse()
        {
            Compute();
            if (mIsSingular) throw new CSException("Can't compute inverse, because matrix is singular");
            int order = this.Rows;
            FloatMatrix inverse = new FloatMatrix(order, order);
            for (int i = 0; i < order; i++)
            {
                inverse.Data[i + order * i] = 1.0f;
            }

            Solve(mMatrix, mPivots, inverse);
            return inverse;
        }

        public override void Clear()
        {
            Array.Clear(mData, 0, mData.Length);
        }

        public override double[] GetRow(int index)
        {
            double[] array = new double[this.Columns];
            for (int i = 0; i < this.Columns; i++) array[i] = ValueAt(index, i);
            return array;
        }

        public override double[] GetColumn(int index)
        {
            double[] array = new double[this.Rows];
            for (int i = 0; i < this.Rows; i++) array[i] = ValueAt(i, index);
            return array;
        }

        public override void SetColumn(int index, double[] source)
        {
            if (index < 0 || index >= Columns) throw new ArgumentOutOfRangeException("index");
            if (source == null) throw new ArgumentNullException("source");
            if (source.Length != Rows) throw new Exception("Array element number is not equal with matrix row count");
            for (int i = 0; i < Rows; i++) ValueAt(i, index, source[i]);
        }

        public override void SetRow(int index, double[] source)
        {
            if (index < 0 || index >= Rows) throw new ArgumentOutOfRangeException("index");
            if (source == null) throw new ArgumentNullException("source");
            if (source.Length != Columns) throw new Exception("Array element number is not equal with matrix Column count");
            for (int j = 0; j < Columns; j++) ValueAt(index, j, source[j]);
        }

        public override Bitmap GetPreview()
        {
            Bitmap ret;
            int i, j;
            int ncols, nrows;
            Int24 pixval = new Int24();
            if (this.Rows * this.Columns > 600 * 600)
            {
                double zoom = Math.Min((double)this.Columns / 600, (double)this.Rows / 600);
                ncols = (int)(zoom * this.Columns);
                nrows = (int)(zoom * this.Rows);
                ret = new Bitmap(ncols, nrows);
                double xfact = (double)this.Columns / ncols;
                double yfact = (double)this.Rows / nrows;

                for (i = 0; i < nrows; i++)
                {
                    for (j = 0; j < ncols; j++)
                    {
                        pixval = (int)this[(int)(Math.Floor(i * yfact)), (int)(Math.Floor(j * xfact))];
                        ret.SetPixel(j, i, pixval.ToColor());
                    }
                }
                return ret;
            }
            else
            {
                ret = new Bitmap(this.Columns, this.Rows);
                for (i = 0; i < this.Rows; i++)
                {
                    for (j = 0; j < this.Columns; j++)
                    {
                        pixval = (int)this[i, j];
                        ret.SetPixel(j, i, pixval.ToColor());
                    }
                }
                return ret;
            }
        }

        #endregion

        #region Features

        public static FloatMatrix operator +(FloatMatrix left, FloatMatrix right)
        {
            if (left == null) throw new ArgumentNullException("left");
            if (right == null) throw new ArgumentNullException("right");
            if (left.Rows != right.Rows || left.Columns != right.Columns) throw new Exception("Input matrix row/column count mismatch");
            FloatMatrix ret = new FloatMatrix(left.Rows, left.Columns);
            for (int i = 0; i < left.Rows; i++)
            {
                for (int j = 0; j < right.Columns; j++)
                {
                    ret[i, j] = left[i, j] + right[i, j];
                }
            }
            return ret;
        }

        public static FloatMatrix operator -(FloatMatrix left, FloatMatrix right)
        {
            if (left == null) throw new ArgumentNullException("left");
            if (right == null) throw new ArgumentNullException("right");
            if (left.Rows != right.Rows || left.Columns != right.Columns) throw new Exception("Input matrix row/column count mismatch");
            FloatMatrix ret = new FloatMatrix(left.Rows, left.Columns);
            for (int i = 0; i < left.Rows; i++)
            {
                for (int j = 0; j < right.Columns; j++)
                {
                    ret[i, j] = left[i, j] - right[i, j];
                }
            }
            return ret;
        }

        public static FloatMatrix operator *(FloatMatrix left, FloatMatrix right)
        {
            if (left == null) throw new ArgumentNullException("left");
            if (right == null) throw new ArgumentNullException("right");
            if (left.Columns != right.Rows) throw new Exception("Input matrix row/column count mismatch");
            FloatMatrix ret = new FloatMatrix(left.Rows, right.Columns);

            for (int j = 0; j != right.Columns; j++)
            {
                for (int i = 0; i != left.Rows; i++)
                {
                    double s = 0;
                    for (int l = 0; l != left.Columns; l++)
                    {
                        s += left[i, l] * right[l, j];
                    }
                    ret[i, j] = s;
                }
            }
            return ret;
        }

        public static FloatMatrix operator /(FloatMatrix left, FloatMatrix right)
        {
            if (left == null) throw new ArgumentNullException("left");
            if (right == null) throw new ArgumentNullException("right");
            if (left.Rows != right.Rows && left.Columns != right.Columns) throw new Exception("Input matrix row/column count mismatch");
            FloatMatrix ret = new FloatMatrix(left.Rows, left.Columns);

            for (int i = 0; i < ret.Rows; i++)
            {
                for (int j = 0; j < ret.Columns; j++)
                {
                    ret[i, j] = left[i, j] / right[i, j];
                }
            }
            return ret;
        }

        #endregion

        protected FloatMatrix mMatrix;
        private bool mComputed;
        private double mDeterminant = Double.MinValue;
        private bool mIsSingular;
        protected int[] mPivots;

        private void Compute()
        {
            if (!mComputed)
            {
                mMatrix = (FloatMatrix)this.Clone();
                mPivots = new int[mMatrix.Rows];

                DoCompute(mMatrix, mPivots);
                for (int j = 0; j < mMatrix.Rows; j++)
                {
                    if (FloatingPointUtilities.EqualsWithinSignificantFigures(mMatrix.ValueAt(j, j), 0.0, 15))
                    {
                        mIsSingular = true;
                        break;
                    }
                }
                mComputed = true;
            }
        }


        private static void Pivot(int m, int n, float[] B, int[] pivots)
        {
            for (int i = 0; i < pivots.Length; i++)
            {
                if (pivots[i] != i)
                {
                    int p = pivots[i];
                    for (int j = 0; j < n; j++)
                    {
                        int indexk = j * m;
                        int indexkp = indexk + p;
                        int indexkj = indexk + i;
                        double temp = B[indexkp];
                        B[indexkp] = B[indexkj];
                        B[indexkj] = (float)temp;
                    }
                }
            }
        }

        private static void Solve(int order, int columns, float[] factor, int[] pivots, float[] data)
        {
            Pivot(order, columns, data, pivots);

            // Solve L*Y = B(piv,:)
            for (int k = 0; k < order; k++)
            {
                int korder = k * order;
                for (int i = k + 1; i < order; i++)
                {
                    for (int j = 0; j < columns; j++)
                    {
                        int index = j * order;
                        data[i + index] -= data[k + index] * factor[i + korder];
                    }
                }
            }
            // Solve U*X = Y;
            for (int k = order - 1; k >= 0; k--)
            {
                int korder = k + k * order;
                for (int j = 0; j < columns; j++)
                {
                    data[k + j * order] /= factor[korder];
                }
                korder = k * order;
                for (int i = 0; i < k; i++)
                {
                    for (int j = 0; j < columns; j++)
                    {
                        int index = j * order;
                        data[i + index] -= data[k + index] * factor[i + korder];
                    }
                }
            }
        }

        protected void Solve(AMatrix factor, int[] pivots, AMatrix result)
        {
            Solve(factor.Rows, result.Columns, ((FloatMatrix)factor).Data, pivots, ((FloatMatrix)result).Data);
        }

        protected void DoCompute(FloatMatrix data, int[] pivots)
        {

            FloatMatrix matrix = data;
            int order = matrix.Rows;
            for (int i = 0; i < order; i++) pivots[i] = i;
            double[] LUcolj = new double[order];
            int indexj, indexjj;

            // Outer loop.
            for (int j = 0; j < order; j++)
            {
                indexj = j * order;
                indexjj = indexj + j;
                // Make a copy of the j-th column to localize references.
                for (int i = 0; i < order; i++) LUcolj[i] = matrix.Data[indexj + i];

                // Apply previous transformations.
                for (int i = 0; i < order; i++)
                {
                    // Most of the time is spent in the following dot product.
                    int kmax = System.Math.Min(i, j);
                    double s = 0.0;
                    for (int k = 0; k < kmax; k++) s += matrix.Data[k * order + i] * LUcolj[k];
                    matrix.Data[indexj + i] = (float)(LUcolj[i] -= s);
                }

                // Find pivot and exchange if necessary.
                int p = j;
                for (int i = j + 1; i < order; i++)
                {
                    if (Math.Abs(LUcolj[i]) > Math.Abs(LUcolj[p])) p = i;
                }
                if (p != j)
                {
                    for (int k = 0; k < order; k++)
                    {
                        int indexk = k * order;
                        int indexkp = indexk + p;
                        int indexkj = indexk + j;
                        double temp = matrix.Data[indexkp];
                        matrix.Data[indexkp] = matrix.Data[indexkj];
                        matrix.Data[indexkj] = (float)temp;
                    }
                    pivots[j] = p;
                }

                // Compute multipliers.
                if (j < order & matrix.Data[indexjj] != 0.0)
                {
                    for (int i = j + 1; i < order; i++) matrix.Data[indexj + i] /= matrix.Data[indexjj];
                }
            }
        }
    }
}
